###############################################################################
#### Indigenous Identity, Natural Resources and Social Conflict in Bolivia ####
####                            Replication:Figures                        ####
###############################################################################

rm(list=ls(all=TRUE))

# Set working directory
setwd(“~/Dropbox/GIGA/ERC/AnalysisBolivia/Replication_Archive")

# load packages
library(raster)
library(maptools)
library(rgdal)
library(raster)
library(rgeos)
library(lattice)
library(fields)
library(gplots)
library(spdep)
library(cshapes)
library(countrycode)
library(sp)
library(Zelig)
library(foreign)
library(ggplot2)
library(scales)
library(plyr)

############################
## Load Data 
############################
bolivia <- read.dta("SocialConflictBolivia.dta")




#############################
## MAKE PLOTS FOR PAPER
#############################

# Figure 1
events <- ddply(bolivia,.(Year_x),summarize,total=sum(freq))

p <- qplot(year, total, data=events, geom="line",xlab="Year",ylab="Conflict Events")
labels <- seq(2000, 2011,1)
labels2 <- seq(0,1000,50)
p <- p + theme_bw()+ scale_x_continuous(breaks=labels, labels=as.character(labels))+ scale_y_continuous(breaks =labels2)
p




# Figure 2
llCRS<-CRS("+proj=longlat +ellps=WGS84")
shape <- readShapePoly("/Users/Jan/Dropbox/GIGA/ERC/AnalysisBolivia/Replication_Archive/bolivia_-_second_administrative_level_boundaries_un_-_salb.shp",proj4string=llCRS)
plot(shape)

events <- ddply(bolivia,~Province,summarise,totevents=sum(freq))
shape@data$Province <- as.factor(c("AITU","ABUN","ALIB","ANIB","ANSA","ANAR","ANQU","ARAN",
                                   "AROM","ARQU","ATAH","EDAV","JMAV","AYOP","JAPA",
                                   "BASA","BELB","BOLI","ECAM","NACA","CAPI","CAVI",
                                   "CARA","CARR","CERC","CADO","CEDO","CECA","CHAP",
                                   "CHAR","CHAY","CHOS","CORD","COSA","DACA","ENBA",
                                   "EARC","FLOR","FRTA","BEBI","FERO","JMPA","JBAL",
                                   "GEBU","GJOR","GRCH","GVIL","GUAR","HESI","ICHI",
                                   "INGA","INQU","ITEN","JMLI","LACA","LARJ","LITO",
                                   "JRLO","LAND","LUCA","MADI","MAMO","MAKA","MMCA",
                                   "MAPI","MARB","PUME","EUME","MIZQ","MOOM","MOXO",
                                   "MUNE","MURI","NISU","NCAR","NCHI","NCIN","NLIP",
                                   "NYUN","NUCH","BUOC","OBSA","OMAS","OROP","PACA",
                                   "PADA","POOP","PUNA","QUIL","RABU","SAJA","SPTO",
                                   "SARA","SAUC","SEPA","SCAR","SCHI","SCIN","SLIP",
                                   "SYUN","TAPA","TIRA","TBAR","TOFR","TOMI","VACD",
                                   "VAGR","JMVE","IGWA","YACU","YAMP","JZUD"))

shape@data <- data.frame(shape@data, events[match(shape@data[,"Province"], events[,"Province"]),])



plot.heat <- function(counties.map,z,title=NULL,breaks=NULL,reverse=FALSE,cex.legend=1,bw=.2,col.vec=NULL,plot.legend=TRUE) {
  ##Break down the value variable
  if (is.null(breaks)) {
    breaks=
      seq(
        floor(min(counties.map@data[,z],na.rm=TRUE)*10)/10
        ,
        ceiling(max(counties.map@data[,z],na.rm=TRUE)*10)/10
        ,.1)
  }
  counties.map@data$zCat <- cut(counties.map@data[,z],breaks,include.lowest=TRUE)
  cutpoints <- levels(counties.map@data$zCat)
  if (is.null(col.vec)) col.vec <- gray.colors(length(levels(counties.map@data$zCat)))
  if (reverse) {
    cutpointsColors <- rev(col.vec)
  } else {
    cutpointsColors <- col.vec
  }
  levels(counties.map@data$zCat) <- cutpointsColors
  plot(counties.map,border="transparent",lwd=0.00001,axes = FALSE, las = 1,col=as.character(counties.map@data$zCat))
  if (plot.legend) legend("bottomleft", cutpoints, fill = cutpointsColors,bty="n",title=title,cex=cex.legend)
  
}

dist <- quantile(shape@data$totevents)
plot.heat(shape,z="totevents",breaks=dist,reverse=TRUE,plot.legend=T)



#######################################
## Figure 3
#######################################
#Substantive Effect, Indigenous Share

bolivia$gasindig <- bolivia$indig_dum*bolivia$gas_count

baseline <- glm.nb(freq~lpop+socecon+lit_dum+mnt2000+capdist
                   +roads_count+gas_count+indig_dum+gasindig+gd2+
                     gd3+gd4+gd5+gd6+gd7+gd8+gd9+gd10+gd11,
                   data=bolivia)
summary(baseline)

means <- coef(baseline)
var <- vcov(baseline)
draws <- mvrnorm(n=1000, means, var)
concvec <- seq(min(bolivia$indig_dum,na.rm=T), max(bolivia$indig_dum,na.rm=T), by=((max(bolivia$indig_dum,na.rm=T)-min(bolivia$indig_dum,na.rm=T))/(500-1)))

data <- bolivia

meanvec <- lowervec <- highervec <- NULL
for(i in 1:length(concvec)){
  out <- draws[,1] + mean(data$lpop,na.rm=T)*draws[,2]+mean(data$socecon,na.rm=T)*draws[,3]+
    mean(data$lit_dum,na.rm=T)*draws[,4]+mean(data$mnt2000,na.rm=T)*draws[,5]+
    mean(data$capdist,na.rm=T)*draws[,6]+mean(data$roads_count,na.rm=T)*draws[,7]+
    0*draws[,8]+concvec[i]*draws[,9]+0*concvec[i]*draws[,10]+1*draws[,16]
  prob <- exp(out)
  
  meanhelp <- means[1] + mean(data$lpop,na.rm=T)*means[2]+mean(data$socecon,na.rm=T)*means[3]+
    mean(data$lit_dum,na.rm=T)*means[4]+mean(data$mnt2000,na.rm=T)*means[5]+
    mean(data$capdist,na.rm=T)*means[6]+mean(data$roads_count,na.rm=T)*means[7]+
    0*means[8]+concvec[i]*means[9]+0*concvec[i]*means[10]+1*means[16]
  
  mean <- exp(meanhelp)
  lower <- quantile(prob,0.05)
  higher <- quantile(prob,0.95)
  meanvec <- c(meanvec,mean)
  lowervec <- c(lowervec,lower)
  highervec <- c(highervec,higher)
}

## with gas

meanvec2 <- lowervec2 <- highervec2 <- NULL
for(i in 1:length(concvec)){
  out <- draws[,1] + mean(data$lpop,na.rm=T)*draws[,2]+mean(data$socecon,na.rm=T)*draws[,3]+
    mean(data$lit_dum,na.rm=T)*draws[,4]+mean(data$mnt2000,na.rm=T)*draws[,5]+
    mean(data$capdist,na.rm=T)*draws[,6]+mean(data$roads_count,na.rm=T)*draws[,7]+
    2*draws[,8]+concvec[i]*draws[,9]+2*concvec[i]*draws[,10]+1*draws[,16]
  prob2 <- exp(out)
  
  meanhelp2 <- means[1] + mean(data$lpop,na.rm=T)*means[2]+mean(data$socecon,na.rm=T)*means[3]+
    mean(data$lit_dum,na.rm=T)*means[4]+mean(data$mnt2000,na.rm=T)*means[5]+
    mean(data$capdist,na.rm=T)*means[6]+mean(data$roads_count,na.rm=T)*means[7]+
    2*means[8]+concvec[i]*means[9]+2*concvec[i]*means[10]+1*means[16]
  
  mean2 <- exp(meanhelp2)
  lower2 <- quantile(prob2,0.05)
  higher2 <- quantile(prob2,0.95)
  meanvec2 <- c(meanvec2,mean2)
  lowervec2 <- c(lowervec2,lower2)
  highervec2 <- c(highervec2,higher2)
}

pdf("IndigenousEffect.pdf",height=7, width=7*1.62)
plot(concvec,meanvec,type="n",ylim=c(0,10), xlab="Share Indigenous Population", ylab="Expected Conflict Events")
polygon(c(concvec, rev(concvec)), c(lowervec,rev(highervec)), col="grey90", border=NA)
points(concvec,meanvec,type="l",lty=1)
polygon(c(concvec, rev(concvec)), c(lowervec2,rev(highervec2)), col="grey65", border=NA)
points(concvec,meanvec2,type="l",lty=1)
rug(data$indig_dum)
dev.off()


#######################################
## Figure X1
#######################################

#Substantive Effect, Gas

baseline <- glm.nb(freq~lpop+socecon+lit_dum+mnt2000+capdist
                   +roads_count+gas_count+indig_dum+gasindig+gd2+
                     gd3+gd4+gd5+gd6+gd7+gd8+gd9+gd10+gd11,
                   data=bolivia)
summary(baseline)

means <- coef(baseline)
var <- vcov(baseline)
draws <- mvrnorm(n=1000, means, var)
concvec <- seq(min(bolivia$gas_count,na.rm=T), max(bolivia$gas_count,na.rm=T), by=1)

data <- bolivia

meanvec <- lowervec <- highervec <- NULL
for(i in 1:length(concvec)){
  out <- draws[,1] + mean(data$lpop,na.rm=T)*draws[,2]+mean(data$socecon,na.rm=T)*draws[,3]+
    mean(data$lit_dum,na.rm=T)*draws[,4]+mean(data$mnt2000,na.rm=T)*draws[,5]+
    mean(data$capdist,na.rm=T)*draws[,6]+mean(data$roads_count,na.rm=T)*draws[,7]+
    concvec[i]*draws[,8]+0.48*draws[,9]+0.48*concvec[i]*draws[,10]+1*draws[,16]
  prob <- exp(out)
  
  meanhelp <- means[1] + mean(data$lpop,na.rm=T)*means[2]+mean(data$socecon,na.rm=T)*means[3]+
    mean(data$lit_dum,na.rm=T)*means[4]+mean(data$mnt2000,na.rm=T)*means[5]+
    mean(data$capdist,na.rm=T)*means[6]+mean(data$roads_count,na.rm=T)*means[7]+
    concvec[i]*means[8]+0.48*means[9]+0.48*concvec[i]*means[10]+1*means[16]
  
  mean <- exp(meanhelp)
  lower <- quantile(prob,0.05)
  higher <- quantile(prob,0.95)
  meanvec <- c(meanvec,mean)
  lowervec <- c(lowervec,lower)
  highervec <- c(highervec,higher)
}

meanvec2 <- lowervec2 <- highervec2 <- NULL
for(i in 1:length(concvec)){
  out <- draws[,1] + mean(data$lpop,na.rm=T)*draws[,2]+mean(data$socecon,na.rm=T)*draws[,3]+
    mean(data$lit_dum,na.rm=T)*draws[,4]+mean(data$mnt2000,na.rm=T)*draws[,5]+
    mean(data$capdist,na.rm=T)*draws[,6]+mean(data$roads_count,na.rm=T)*draws[,7]+
    concvec[i]*draws[,8]+0.98*draws[,9]+0.98*concvec[i]*draws[,10]+1*draws[,16]
  prob2 <- exp(out)
  
  meanhelp2 <- means[1] + mean(data$lpop,na.rm=T)*means[2]+mean(data$socecon,na.rm=T)*means[3]+
    mean(data$lit_dum,na.rm=T)*means[4]+mean(data$mnt2000,na.rm=T)*means[5]+
    mean(data$capdist,na.rm=T)*means[6]+mean(data$roads_count,na.rm=T)*means[7]+
    concvec[i]*means[8]+0.98*means[9]+0.98*concvec[i]*means[10]+1*means[16]
  
  mean2 <- exp(meanhelp2)
  lower2 <- quantile(prob2,0.05)
  higher2 <- quantile(prob2,0.95)
  meanvec2 <- c(meanvec2,mean2)
  lowervec2 <- c(lowervec2,lower2)
  highervec2 <- c(highervec2,higher2)
}



pdf("GasEffect.pdf",height=7, width=7*1.62)
plot(concvec,meanvec,type="n",ylim=c(0,10), xlab="Gas", ylab="Expected Conflict Events")
polygon(c(concvec, rev(concvec)), c(lowervec,rev(highervec)), col="grey90", border=NA)
points(concvec,meanvec,type="l",lty=1)
polygon(c(concvec, rev(concvec)), c(lowervec2,rev(highervec2)), col="grey65", border=NA)
points(concvec,meanvec2,type="l",lty=1)
rug(data$gas_count)
dev.off()


